Basic scMultiome data analysis

Vignette for integrated Seurat/Signac analysis: https://stuartlab.org/signac/articles/pbmc_multiomic.html

Overview

  1. Seurat objects
    a) Create RNA assays
    b) Create ATAC assays
    c) Cell statistics
  2. Filtering
    a) Cell filtering by HQ cells from single omics
    b) Sample statistics
    c) Cell statistics
  3. Preprocessing
    a) ATAC
    b) RNA
  4. Dimensionality reduction and single-cell embedding
    a) ATAC
    b) RNA
    c) Neighbor graph
    d) Non-linear dimension reduction
  5. Save results

Input data

  1. 10x Genomics Multiome-scRNA/ATAC-seq of PBMCs
    PBMC_Multiome_10xTn5_CellRangerARC_fragments.tsv.gz
    PBMC_Multiome_10xTn5_CellRangerARC/raw_feature_bc_matrix
  2. 10x Genomics Multiome-scRNA/TurboATAC-seq (1:6 concentration of Tn5) of PBMCs
    PBMC_Multiome_Tn5hc_50perc_CellRangerARC_fragments.tsv.gz
    PBMC_Multiome_Tn5hc_50perc_CellRangerARC/raw_feature_bc_matrix
  3. 10x Genomics Multiome-scRNA/TurboATAC-seq (1:3 concentration of Tn5) of PBMCs
    PBMC_Multiome_Tn5hc_CellRangerARC_fragments.tsv.gz
    PBMC_Multiome_Tn5hc_CellRangerARC/raw_feature_bc_matrix
  4. HQ cells from ATAC
    RawData_scATAC_*.csv
    For prep see 20230912_PBMC_Multiome_scATAC_UMAP.ipynb
  5. HQ cells from RNA
    RawData_scRNA_*.csv
    For prep see 20230912_PBMC_Multiome_scRNA_UMAP.ipynb

Notebook

20230912_PBMC_Multiome_Co-embedding.ipynb


0. Load libraries, scripts and data

### Libraries
library(magrittr)
library(qs)
library(Signac)
library(Seurat)
library(EnsDb.Hsapiens.v86)
library(BSgenome.Hsapiens.UCSC.hg38)
library(ggplot2)
library(ggpubr)
### Set seed
set.seed(13)
library(future)
plan("multicore", workers = 20)
options(future.globals.maxSize = 500000 * 1024^2) # for 500 Gb RAM
### Samples
samples <- c("PBMC_scATAC", "PBMC_scTurboATAC_16", "PBMC_scTurboATAC_13")

### Sample palette
sample_palette <- c("grey", "blue", "darkblue")
names(sample_palette) <- samples
### Get gene annotations for hg38
annotation <- GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v86)
seqlevelsStyle(annotation) <- "UCSC"
### Data
inputFiles <- c("PBMC_Multiome_10xTn5_CellRangerARC_fragments.tsv.gz",
                "PBMC_Multiome_Tn5hc_50perc_CellRangerARC_fragments.tsv.gz",
                "PBMC_Multiome_Tn5hc_CellRangerARC_fragments.tsv.gz")
names(inputFiles) <- samples
### Filtered ATAC cells
atac_cells <- list()

for (sample in samples){
    atac_cells[[sample]] <- read.table(paste0("RawData_scATAC_", 
                                              sample, ".csv"), sep = ",", header = TRUE)
    colnames(atac_cells[[sample]])[1] <- "CellID_ArchR"
    atac_cells[[sample]]$CellID <- gsub(paste0(sample, "#"), "", atac_cells[[sample]]$CellID_ArchR)
}
### Filtered RNA cells
rna_cells <- list()

for (sample in samples){
    rna_cells[[sample]] <- read.table(paste0("RawData_scRNA_", 
                                              sample, ".csv"), sep = ",", header = TRUE)
    colnames(rna_cells[[sample]])[1] <- "CellID_Seurat"
    rna_cells[[sample]]$CellID <- gsub(paste0(sample, "_"), "", rna_cells[[sample]]$CellID_Seurat)
}

1. Seurat objects

overview

a) Create RNA assays

input_data <- list()
for (sample in samples){
    input_data[[sample]] <- Read10X(data.dir = paste0(sample, "/raw_feature_bc_matrix"))
}
proj_list <- list()
for (sample in samples){
    proj_list[[sample]] <- CreateSeuratObject(counts = input_data[[sample]][["Gene Expression"]], project = sample, assay = "RNA")
    proj_list[[sample]]$orig.ident <- sample
}
for (sample in samples){
    proj_list[[sample]]$percent.mt <- PercentageFeatureSet(proj_list[[sample]], pattern = "^MT-")
}

b) Create ATAC assays

frag_list <- list()

for (sample in samples){
    frag_list[[sample]] <- CreateFragmentObject(inputFiles[sample])
}
atac_list <- list()

for (sample in samples){
    proj_list[[sample]][["ATAC"]] <- CreateChromatinAssay(counts = input_data[[sample]][["Peaks"]], fragments = frag_list[[sample]],
                                                          genome = 'hg38', sep = c(":", "-"), annotation = annotation)
}
for (sample in samples){
    DefaultAssay(proj_list[[sample]]) <- "ATAC"
    
    proj_list[[sample]] <- TSSEnrichment(proj_list[[sample]])
    proj_list[[sample]] <- NucleosomeSignal(proj_list[[sample]])
}
proj_list
$PBMC_scATAC
An object of class Seurat 
156885 features across 720228 samples within 2 assays 
Active assay: ATAC (120284 features, 0 variable features)
 1 other assay present: RNA

$PBMC_scTurboATAC_16
An object of class Seurat 
149742 features across 716191 samples within 2 assays 
Active assay: ATAC (113141 features, 0 variable features)
 1 other assay present: RNA

$PBMC_scTurboATAC_13
An object of class Seurat 
137725 features across 719588 samples within 2 assays 
Active assay: ATAC (101124 features, 0 variable features)
 1 other assay present: RNA

c) Cell statistics

plot_list <- list()

for (sample in samples){
    plot_list[[paste0(sample, "_nCountRNA_nCountATAC")]] <- FeatureScatter(proj_list[[sample]], feature1 = "nCount_RNA", feature2 = "nCount_ATAC", 
                                                                           cols = sample_palette, group.by = "orig.ident") + 
                                                                theme(legend.position="none") + 
                                                                scale_x_log10(labels = scales::comma) + scale_y_log10(labels = scales::comma) +
                                                                ggtitle(paste(sample, " (", length(proj_list[[sample]]$orig.ident), "cells)"))
    plot_list[[paste0(sample, "_percentMT_TSSenrichment")]] <- FeatureScatter(proj_list[[sample]], feature1 = "percent.mt", feature2 = "TSS.enrichment",  
                                                                              cols = sample_palette, group.by = "orig.ident") + 
                                                                theme(legend.position="none") + 
                                                                ggtitle(paste(sample, " (", length(proj_list[[sample]]$orig.ident), "cells)"))
    plot_list[[paste0(sample, "_percentMT_NucleosomeRatio")]] <- FeatureScatter(proj_list[[sample]], feature1 = "percent.mt", feature2 = "nucleosome_signal",  
                                                                              cols = sample_palette, group.by = "orig.ident") + 
                                                                theme(legend.position="none") + 
                                                                ggtitle(paste(sample, " (", length(proj_list[[sample]]$orig.ident), "cells)"))
}
options(repr.plot.width=18, repr.plot.height=5*length(samples))
ggarrange(plotlist=plot_list, ncol=3, nrow=length(samples))
Warning message:
"Transformation introduced infinite values in continuous x-axis"
Warning message:
"Transformation introduced infinite values in continuous y-axis"
Warning message:
"Removed 149665 rows containing missing values (`geom_scattermore()`)."
Warning message:
"Removed 219953 rows containing missing values (`geom_scattermore()`)."
Warning message:
"Transformation introduced infinite values in continuous x-axis"
Warning message:
"Transformation introduced infinite values in continuous y-axis"
Warning message:
"Removed 236481 rows containing missing values (`geom_scattermore()`)."
Warning message:
"Removed 291334 rows containing missing values (`geom_scattermore()`)."
Warning message:
"Transformation introduced infinite values in continuous x-axis"
Warning message:
"Transformation introduced infinite values in continuous y-axis"
Warning message:
"Removed 247437 rows containing missing values (`geom_scattermore()`)."
Warning message:
"Removed 291442 rows containing missing values (`geom_scattermore()`)."
summary <- lapply(samples, 
                       function(sample){c(sample, '01_raw',
                                     nrow(proj_list[[sample]]@meta.data),
                                     median(proj_list[[sample]]$nFeature_RNA, na.rm = TRUE) %>% round(., 0),
                                     median(proj_list[[sample]]$nCount_RNA, na.rm = TRUE) %>% round(., 0),
                                     median(proj_list[[sample]]$percent.mt, na.rm = TRUE) %>% round(., 2),
                                     median(proj_list[[sample]]$nFeature_ATAC, na.rm = TRUE) %>% round(., 0),
                                     median(proj_list[[sample]]$nCount_ATAC, na.rm = TRUE) %>% round(., 0),
                                     median(proj_list[[sample]]$TSS.enrichment, na.rm = TRUE) %>% round(., 2),
                                     median(proj_list[[sample]]$nucleosome_signal, na.rm = TRUE) %>% round(., 2))}) %>% do.call(rbind, .)
colnames(summary) <- c('sample', 'step', 'nCells', 'nFeature_RNA', 'nCount_RNA', 'percent.mt', 'nFeature_ATAC', 'nCount_ATAC', 'TSSenrichment', 'NucleosomeSignal')
summary
A matrix: 3 × 10 of type chr
samplestepnCellsnFeature_RNAnCount_RNApercent.mtnFeature_ATACnCount_ATACTSSenrichmentNucleosomeSignal
PBMC_scATAC 01_raw7202282222.22240.041.29
PBMC_scTurboATAC_1601_raw7161911120.69120 1.75
PBMC_scTurboATAC_1301_raw7195881112.86120 2

2. Filtering

overview
HQ_cells <- data.frame(Sample = samples,
                       HQ_atac_cells = sapply(atac_cells, nrow),
                       HQ_rna_cells = sapply(rna_cells, nrow),
                       HQ_cells = sapply(samples, function(sample){sum(atac_cells[[sample]]$CellID %in% rna_cells[[sample]]$CellID)}))
HQ_cells
A data.frame: 3 × 4
SampleHQ_atac_cellsHQ_rna_cellsHQ_cells
<chr><int><int><int>
PBMC_scATACPBMC_scATAC 414142773620
PBMC_scTurboATAC_16PBMC_scTurboATAC_16678573895999
PBMC_scTurboATAC_13PBMC_scTurboATAC_13576565295077

a) Cell filtering by HQ cells from single omics

for (sample in samples){
    proj_list[[sample]]$HQ_cell <- rownames(proj_list[[sample]]@meta.data)  %in% atac_cells[[sample]]$CellID[which(atac_cells[[sample]]$CellID %in% rna_cells[[sample]]$CellID)]
    proj_list[[sample]] <- subset(proj_list[[sample]], subset = HQ_cell)
}
summary <- lapply(samples, 
                       function(sample){c(sample, '02_HQcells',
                                     nrow(proj_list[[sample]]@meta.data),
                                     median(proj_list[[sample]]$nFeature_RNA) %>% round(., 0),
                                     median(proj_list[[sample]]$nCount_RNA) %>% round(., 0),
                                     median(proj_list[[sample]]$percent.mt) %>% round(., 2),
                                     median(proj_list[[sample]]$nFeature_ATAC) %>% round(., 0),
                                     median(proj_list[[sample]]$nCount_ATAC) %>% round(., 0),
                                     median(proj_list[[sample]]$TSS.enrichment) %>% round(., 2),
                                     median(proj_list[[sample]]$nucleosome_signal) %>% round(., 2))}) %>% do.call(rbind, .) %>% rbind(summary, .)
colnames(summary) <- c('sample', 'step', 'nCells', 'nFeature_RNA', 'nCount_RNA', 'percent.mt', 'nFeature_ATAC', 'nCount_ATAC', 'TSSenrichment', 'NucleosomeSignal')
summary
A matrix: 6 × 10 of type chr
samplestepnCellsnFeature_RNAnCount_RNApercent.mtnFeature_ATACnCount_ATACTSSenrichmentNucleosomeSignal
PBMC_scATAC 01_raw 7202282 2 22.222 4 0.041.29
PBMC_scTurboATAC_1601_raw 7161911 1 20.691 2 0 1.75
PBMC_scTurboATAC_1301_raw 7195881 1 12.861 2 0 2
PBMC_scATAC 02_HQcells3620 184440247.98 15708439304.451.09
PBMC_scTurboATAC_1602_HQcells5999 129322535.64 11426281554.2 1.18
PBMC_scTurboATAC_1302_HQcells5077 139225014.71 11431271803.621.46

b) Sample statistics

plot_list <- list()

for (feature in c("nFeature_ATAC", "nCount_ATAC", "TSS.enrichment", "nucleosome_signal", "nFeature_RNA", "nCount_RNA", "percent.mt")){
    plot_list[[feature]] <- ggplot(lapply(proj_list, function(x){data.frame(Sample = factor(x$orig.ident, levels = samples), 
                                                                            y = x@meta.data[, feature])}) %>% do.call(rbind, .), 
                                   aes(x = Sample, y = y)) + 
                            geom_violin(aes(fill = Sample), alpha = 0.8) + scale_fill_manual(values = sample_palette) +
                            geom_boxplot(alpha = 0.001) + 
                            xlab("") + ylab(feature) +
                            theme_minimal() + theme(legend.position = "none", plot.title = element_text(size=14), 
                                                    axis.title = element_text(size=12), axis.text = element_text(size=12), 
                                                    axis.text.x=element_text(angle=45, vjust=1, hjust=1), 
                                                    legend.title = element_text(size=12), legend.text = element_text(size=12))
}

plot_list$nCount_ATAC <- plot_list$nCount_ATAC + scale_y_log10(labels = scales::comma) + ylab("log10(nCount_ATAC)")
plot_list$nCount_RNA <- plot_list$nCount_RNA + scale_y_log10(labels = scales::comma) + ylab("log10(nCount_RNA)")
options(repr.plot.width=15, repr.plot.height=12)
ggarrange(plotlist = plot_list, ncol = 4, nrow = 2)

c) Cell statistics

plot_list <- list()

for (sample in samples){
    plot_list[[paste0(sample, "_nCountRNA_nCountATAC")]] <- FeatureScatter(proj_list[[sample]], feature1 = "nCount_RNA", feature2 = "nCount_ATAC", 
                                                                           cols = sample_palette, group.by = "orig.ident") + 
                                                                theme(legend.position="none") + 
                                                                scale_x_log10(labels = scales::comma) + scale_y_log10(labels = scales::comma) +
                                                                ggtitle(paste(sample, " (", length(proj_list[[sample]]$orig.ident), "cells)"))
    plot_list[[paste0(sample, "_percentMT_TSSenrichment")]] <- FeatureScatter(proj_list[[sample]], feature1 = "percent.mt", feature2 = "TSS.enrichment",  
                                                                              cols = sample_palette, group.by = "orig.ident") + 
                                                                theme(legend.position="none") + 
                                                                ggtitle(paste(sample, " (", length(proj_list[[sample]]$orig.ident), "cells)"))
    plot_list[[paste0(sample, "_percentMT_NucleosomeRatio")]] <- FeatureScatter(proj_list[[sample]], feature1 = "percent.mt", feature2 = "nucleosome_signal",  
                                                                              cols = sample_palette, group.by = "orig.ident") + 
                                                                theme(legend.position="none") + 
                                                                ggtitle(paste(sample, " (", length(proj_list[[sample]]$orig.ident), "cells)"))
}
options(repr.plot.width=18, repr.plot.height=5*length(samples))
ggarrange(plotlist=plot_list, ncol=3, nrow=length(samples))

3. Preprocessing

overview

a) ATAC

for (sample in samples){
    DefaultAssay(proj_list[[sample]]) <- "ATAC"
    
    proj_list[[sample]] <- FindTopFeatures(proj_list[[sample]], min.cutoff = 5)
    proj_list[[sample]] <- RunTFIDF(proj_list[[sample]])
}

b) RNA

for (sample in samples){
    DefaultAssay(proj_list[[sample]]) <- "RNA"
    
    proj_list[[sample]] <- SCTransform(proj_list[[sample]])
}

4. Dimensionality reduction and single-cell embedding

overview

a) ATAC

for (sample in samples){
    DefaultAssay(proj_list[[sample]]) <- "ATAC"
    
    proj_list[[sample]] <- RunSVD(proj_list[[sample]])
}
plot_list <- list()

for (sample in samples){
    plot_list[[sample]] <- ElbowPlot(proj_list[[sample]], reduction = "lsi", ndims = 50)
}
options(repr.plot.width=6*length(samples), repr.plot.height=5)
ggarrange(plotlist = plot_list, ncol = length(samples), nrow = 1)
LSIs_list <- list(25, 25, 25)
names(LSIs_list) <- samples

b) RNA

for (sample in samples){
    DefaultAssay(proj_list[[sample]]) <- "SCT"
    
    proj_list[[sample]] <- RunPCA(proj_list[[sample]])
}
plot_list <- list()

for (sample in samples){
    plot_list[[sample]] <- ElbowPlot(proj_list[[sample]], reduction = "pca", ndims = 50)
}
options(repr.plot.width=6*length(samples), repr.plot.height=5)
ggarrange(plotlist = plot_list, ncol = length(samples), nrow = 1)
PCs_list <- list(30, 30, 30)
names(PCs_list) <- samples

c) Neighbor graph

for (sample in samples){
    proj_list[[sample]] <- FindMultiModalNeighbors(proj_list[[sample]], reduction.list = list("pca", "lsi"), 
                                                   dims.list = list(1:PCs_list[[sample]], 2:LSIs_list[[sample]]), 
                                                   modality.weight.name = "RNA.weight")
    
    for (resolution in c(0.2, 0.5, 0.8)){
        proj_list[[sample]] <- FindClusters(proj_list[[sample]], resolution = resolution, graph.name = "wsnn")
    }
}

d) Non-linear dimension reduction

for (sample in samples){
    proj_list[[sample]] <- RunUMAP(proj_list[[sample]], nn.name = "weighted.nn", assay = "RNA")
}
plot_list <- list()

for (sample in samples){
    plot_list[[sample]] <- list()
    
    plot_list[[sample]]$clusters0.2 <- DimPlot(proj_list[[sample]], reduction = "umap", group.by = "wsnn_res.0.2")
    plot_list[[sample]]$clusters0.5 <- DimPlot(proj_list[[sample]], reduction = "umap", group.by = "wsnn_res.0.5")
    plot_list[[sample]]$clusters0.8 <- DimPlot(proj_list[[sample]], reduction = "umap", group.by = "wsnn_res.0.8")
    
    plot_list[[sample]]$nCount_RNA <- FeaturePlot(proj_list[[sample]], reduction = "umap", features = "nCount_RNA")
    plot_list[[sample]]$nFeature_RNA <- FeaturePlot(proj_list[[sample]], reduction = "umap", features = "nFeature_RNA")
    plot_list[[sample]]$percentMT <- FeaturePlot(proj_list[[sample]], reduction = "umap", features = "percent.mt")
    
    plot_list[[sample]]$nCount_ATAC <- FeaturePlot(proj_list[[sample]], reduction = "umap", features = "nCount_ATAC")
    plot_list[[sample]]$nFeature_ATAC <- FeaturePlot(proj_list[[sample]], reduction = "umap", features = "nFeature_ATAC")
    plot_list[[sample]]$NucRatio <- FeaturePlot(proj_list[[sample]], reduction = "umap", features = "nucleosome_signal")
    plot_list[[sample]]$TSSenrichment <- FeaturePlot(proj_list[[sample]], reduction = "umap", features = "TSS.enrichment")

}
options(repr.plot.width=20, repr.plot.height=24)
for (sample in samples){
    ggarrange(plotlist = plot_list[[sample]], ncol = 3, nrow = 4) %>% print(.)
}

5. Save results

### Save Seurat object
qsave(proj_list, 'SeuratObjects.qs')
### Save Sample statistics
write.table(summary, 'Summary.csv', sep=',', col.names = T, row.names = F)
### Save raw data
lapply(samples, function(x){all(rownames(proj_list[[sample]]@meta.data) == rownames(proj_list[[sample]]@reductions$umap@cell.embeddings)) %>% print(.);
                            df <- cbind(proj_list[[sample]]@meta.data,
                                        proj_list[[sample]]@reductions$umap@cell.embeddings)
                            write.csv(df, paste0("RawData_", x, ".csv"))})

#

Isabelle Seufert
Division of Chromatin Networks, DKFZ
12.09.2023

sessionInfo()
R version 4.2.2 (2022-10-31)
Platform: x86_64-conda-linux-gnu (64-bit)
Running under: CentOS Linux 7 (Core)

Matrix products: default
BLAS/LAPACK: /home/bq_ilander/miniconda3/envs/RWireX_v0.2.02/lib/libopenblasp-r0.3.21.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats4    stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] future_1.32.0                     ggpubr_0.4.0                     
 [3] ggplot2_3.4.2                     BSgenome.Hsapiens.UCSC.hg38_1.4.4
 [5] BSgenome_1.66.3                   rtracklayer_1.58.0               
 [7] Biostrings_2.66.0                 XVector_0.38.0                   
 [9] EnsDb.Hsapiens.v86_2.99.0         ensembldb_2.22.0                 
[11] AnnotationFilter_1.22.0           GenomicFeatures_1.50.2           
[13] AnnotationDbi_1.60.0              Biobase_2.58.0                   
[15] GenomicRanges_1.50.1              GenomeInfoDb_1.34.9              
[17] IRanges_2.32.0                    S4Vectors_0.36.0                 
[19] BiocGenerics_0.44.0               SeuratObject_4.1.3               
[21] Seurat_4.3.0                      Signac_1.9.0                     
[23] qs_0.25.4                         magrittr_2.0.3                   

loaded via a namespace (and not attached):
  [1] utf8_1.2.3                  R.utils_2.12.1             
  [3] spatstat.explore_3.1-0      reticulate_1.26            
  [5] tidyselect_1.2.0            RSQLite_2.3.1              
  [7] htmlwidgets_1.6.2           grid_4.2.2                 
  [9] BiocParallel_1.32.5         Rtsne_0.16                 
 [11] munsell_0.5.0               ragg_1.2.5                 
 [13] codetools_0.2-19            ica_1.0-3                  
 [15] pbdZMQ_0.3-9                miniUI_0.1.1.1             
 [17] withr_2.5.0                 spatstat.random_3.1-4      
 [19] colorspace_2.1-0            progressr_0.13.0           
 [21] filelock_1.0.2              knitr_1.42                 
 [23] uuid_1.1-0                  rstudioapi_0.14            
 [25] ROCR_1.0-11                 ggsignif_0.6.4             
 [27] tensor_1.5                  listenv_0.9.0              
 [29] labeling_0.4.2              MatrixGenerics_1.10.0      
 [31] repr_1.1.6                  GenomeInfoDbData_1.2.9     
 [33] polyclip_1.10-4             farver_2.1.1               
 [35] bit64_4.0.5                 parallelly_1.35.0          
 [37] vctrs_0.6.2                 generics_0.1.3             
 [39] xfun_0.39                   biovizBase_1.46.0          
 [41] BiocFileCache_2.6.0         R6_2.5.1                   
 [43] bitops_1.0-7                spatstat.utils_3.0-3       
 [45] cachem_1.0.8                DelayedArray_0.24.0        
 [47] promises_1.2.0.1            BiocIO_1.8.0               
 [49] scales_1.2.1                nnet_7.3-19                
 [51] gtable_0.3.3                Cairo_1.6-0                
 [53] globals_0.16.2              goftest_1.2-3              
 [55] rlang_1.1.1                 systemfonts_1.0.4          
 [57] RcppRoll_0.3.0              splines_4.2.2              
 [59] rstatix_0.7.2               lazyeval_0.2.2             
 [61] dichromat_2.0-0.1           checkmate_2.2.0            
 [63] broom_1.0.4                 spatstat.geom_3.2-1        
 [65] yaml_2.3.7                  reshape2_1.4.4             
 [67] abind_1.4-5                 backports_1.4.1            
 [69] httpuv_1.6.10               Hmisc_5.1-0                
 [71] tools_4.2.2                 ellipsis_0.3.2             
 [73] RColorBrewer_1.1-3          ggridges_0.5.4             
 [75] Rcpp_1.0.10                 plyr_1.8.8                 
 [77] base64enc_0.1-3             progress_1.2.2             
 [79] zlibbioc_1.44.0             purrr_1.0.1                
 [81] RCurl_1.98-1.9              prettyunits_1.1.1          
 [83] rpart_4.1.19                deldir_1.0-6               
 [85] pbapply_1.7-0               cowplot_1.1.1              
 [87] zoo_1.8-12                  SummarizedExperiment_1.28.0
 [89] ggrepel_0.9.3               cluster_2.1.4              
 [91] data.table_1.14.8           scattermore_1.0            
 [93] lmtest_0.9-40               RANN_2.6.1                 
 [95] ProtGenerics_1.30.0         fitdistrplus_1.1-11        
 [97] matrixStats_0.63.0          stringfish_0.15.7          
 [99] hms_1.1.3                   patchwork_1.1.2            
[101] mime_0.12                   evaluate_0.21              
[103] xtable_1.8-4                XML_3.99-0.14              
[105] gridExtra_2.3               compiler_4.2.2             
[107] biomaRt_2.54.0              tibble_3.2.1               
[109] KernSmooth_2.23-21          crayon_1.5.2               
[111] R.oo_1.25.0                 htmltools_0.5.5            
[113] later_1.3.1                 Formula_1.2-5              
[115] tidyr_1.3.0                 RcppParallel_5.1.6         
[117] DBI_1.1.3                   RApiSerialize_0.1.2        
[119] dbplyr_2.3.2                MASS_7.3-60                
[121] rappdirs_0.3.3              car_3.1-2                  
[123] Matrix_1.5-4                cli_3.6.1                  
[125] R.methodsS3_1.8.2           parallel_4.2.2             
[127] igraph_1.4.2                pkgconfig_2.0.3            
[129] GenomicAlignments_1.34.0    foreign_0.8-84             
[131] sp_1.6-0                    IRdisplay_1.1              
[133] plotly_4.10.1               spatstat.sparse_3.0-1      
[135] xml2_1.3.3                  VariantAnnotation_1.44.0   
[137] stringr_1.5.0               digest_0.6.31              
[139] sctransform_0.3.5           RcppAnnoy_0.0.20           
[141] spatstat.data_3.0-1         rmarkdown_2.21             
[143] leiden_0.4.3                fastmatch_1.1-3            
[145] htmlTable_2.4.1             uwot_0.1.14                
[147] restfulr_0.0.15             curl_4.3.3                 
[149] shiny_1.7.4                 Rsamtools_2.14.0           
[151] rjson_0.2.21                lifecycle_1.0.3            
[153] nlme_3.1-162                jsonlite_1.8.4             
[155] carData_3.0-5               viridisLite_0.4.2          
[157] fansi_1.0.4                 pillar_1.9.0               
[159] lattice_0.21-8              KEGGREST_1.38.0            
[161] fastmap_1.1.1               httr_1.4.6                 
[163] survival_3.5-5              glue_1.6.2                 
[165] png_0.1-8                   bit_4.0.5                  
[167] stringi_1.7.12              blob_1.2.4                 
[169] textshaping_0.3.6           memoise_2.0.1              
[171] IRkernel_1.3.1              dplyr_1.1.2                
[173] irlba_2.3.5.1               future.apply_1.10.0